Ginzburg-Landau Vortex Lattice in Superconductor Films of Finite Thickness 
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The Ginzburg-Landau equations are solved for ideally periodic vortex lattices in superconducting 
films of arbitrary thickness in a perpendicular magnetic field. The order parameter, current density, 
magnetic moment, and the 3-dimensional magnetic field inside and outside the film are obtained in 
the entire ranges of the applied magnetic field, Ginzburg Landau parameter re, and film thickness. 
The superconducting order parameter varies very little near the surface (~ 1%) and the energy 
of the film surface is small. The shear modulus c^e of the triangular vortex lattice in thin films 
coincides with the bulk C66 taken at large re. In thin type-I superconductor films with re < 1/ y/2, cee 
can be positive at low fields and negative at high fields. The magnetization of thin films at small 
applied fields is enhanced beyond its bulk value —H c \ due to the energy of the magnetic stray field. 

PACS numbers: 74.25.Qt, 74.78.Bz, 74.78.Db 



I. INTRODUCTION 

Since Abrikosov'si prediction of the flux-line lattice 
in Type-II superconductors from Ginzburg-Landau (GL) 
theory, most theoretical work on this vortex lattice in 
bulk and thin film superconductors considered the sit- 
uation when the applied magnetic field and the aver- 
age induction B are close to the upper critical field 
B C 2 = /J-qHc2i since analytical solutions may be obtained 
for this particular case. In the opposite limit of low in- 
duction B <§; B C 2, the properties of an isolated vortex 
and the interaction between vortices are obtained to good 
approximation from the London theory when the GL pa- 
rameter k is not too small, 2re 2 ^> The problem 
of an isolated vortex in thin films was solved analyti- 
cally within London theory by Pearl^; the interaction 
energy of such Pearl vortices (or pancake vortices^) is 
easily calculated by noting that within London theory 
the currents and magnetic fields of the vortices superim- 
pose linearly and that the force on a vortex equals the 
thickness-integrated super-current density at the vortex 
core times the quantum of flux $o- I n thin films with 
thickness d smaller than the London magnetic penetra- 
tion depth A the range of the vortex-vortex interaction 
is increased to the effective penetration depth A = X 2 /d 
since the interaction now occurs mainly via the magnetic 
stray field outside the film^ Vortices in superconduct- 
ing films of finite thickness [d < A and d > A) and in the 
superconducting half space (d 3> A) were calculated from 
GL theory^ and London theoryi^^ 

At larger reduced induction b = B/B C 2 > 0.05 when 
the London theory does not apply, the properties of the 
GL vortex lattice have to be computed numerically. A 
very efficient method^ uses Fourier series as trial func- 
tions for the GL function |-0(a;, y)| 2 and magnetic field 
B(x, y) and minimizes the GL free energy with respect 
to a finite number of Fourier coefficients. This numeri- 
cal method was recently improvecUii^ by solving the GL 
equations iteratively with high precision. 

The present paper extends this two-dimensional (2D) 
method to the 3D problem of a film of arbitrary thickness 



containing a periodic lattice of GL vortices oriented per- 
pendicular to the film plane. Due to the Fourier ansatz, 
the magnetic stray field energy is easily accounted for 
in this method. Moreover, it turns out that the exten- 
sion from 2D to a 3D problem only slightly increases 
the required computation time and computer memory, 
so that high precision can be achieved easily on a Per- 
sonal Computer. Like in Refs. [T1IIT2I we consider here 
vortex lattices with arbitrary shape of the unit cell con- 
taining one vortex, i.e., our method computes triangu- 
lar, rectangular, square lattices, etc., and yields also the 
shear moduli^ of the equilibrium lattices. The approx- 
imate shear modulus cqq of the triangular vortex lattice 
in thin films was com put ed from GL theory for b <C 1 
and 1 — b <C 1 in Ref. [lj. For early work on films with 
perpendicular vortex lattice see Refs. I2l3lll5lll6lll7lll8lll9l 

Though we consider here isotropic superconductors, 
the corresponding results for anisotropic superconduc- 
tors with principal symmetry axes along x,y,z may be 
obtained from this isotropic method by scaling the co- 
ordinates and introducing an effective GL parameter 
- 20,21,22 -p nc ma g ne tic field of a vortex inside a uniax- 
ially anisotropic superconductor with surface parallel to 
the a, c symmetry plane and perpendicular to the vortex 
line was calculated from anisotropic London theory^ and 
compared with experiments in Ref. 12.1 

The main effect of the flat surface of a superconductor 
film or half space is the widening of the magnetic field 
lines of the vortices as they approach the surface. This 
widening minimizes the sum of the bulk free energy plus 
the energy of the magnetic stray field outside the super- 
conductor. The resulting magnetic field lines cross the 
superconductor surface smoothly^ee Fig. 1 for the vor- 
tex lattice and Figs. 1,2 of Ref. [3 for isolated vortices. 
One can see that for the vortex lattice the field lines at 
the boundary of the Wigner-Seitz cell are exactly par- 
allel to z, inside and outside the superconductor, and at 
some distance outside from the surface (*=s half the vortex 
spacing) the magnetic field becomes uniform and thus the 
field lines are parallel and equidistant. For the isolated 
vortex, the field lines away from the surface approach ra- 
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dial lines as if they would originate from a point source, 
a magnetic monopole with strength 2$ situated on the 
vortex core at a distance 1.27A below the surfaced 

In Ref . |a the widening of the field lines inside the su- 
perconductor was missed, but some modification of the 
superconductor order parameter near the surface was cal- 
culated from GL theory. Below we obtain that the correct 
modification of \ip\ 2 is very small: the vortex core, visu- 
alized as contour lines of \tp(x, y, z)\ 2 , widens near the 
surface by only a few percent. 

The outline of this paper is as follows. In Set. 2 the 
solution method is outlined. Section 3 presents a selec- 
tion of results for thin and thick films: Magnetic field 
lines, profiles of the order parameter and magnetic field, 
the variances of the periodic order parameter and mag- 
netic field inside and outside the film, surface energy and 
stray-field energy, and shear modulus of the triangular 
vortex lattice in the film. A summary is given in Set. 4. 



II. SOLUTION METHOD 

The properties of the vortex lattice within GL theory 
are calculated by minimizing the GL free energy of the 
superconductor with respect to the complex GL function 
ip(r) and to the vector potential A(r) of the local mag- 
netic induction B(r) = V x A. In the usual reduced 
unitsi^ (length A, magnetic induction v2~Mo-^o energy 
density /io-ff 2 , where H c = i7 c2 /(v / 2^) is the thermo- 
dynamic critical field) the spatially averaged free energy 
density F of the GL theory referred to the Meissner state 

= 1, B = 0) within the superconductor reads 



F = 



2\2 



-aU 



+ B 2 



(1) 



Here (...) = (1/V) J v d 3 r . . . means spatial averaging 
over the superconductor with volume V. Introducing 
the supervelocity Q(r) — A — S7tp/n and the magnitude 
/(r) = \ip\ of ^i(r) = /(r) exp[«/?(r)] one may write F as 
a functional of the real and gauge-invariant functions / 
and Q, 



F 



(I-/ 2 ) 2 , (V/) 2 



/ 2 Q 2 + (VxQ) ; 



(2) 



In the presence of vortices Q(r) has to be chosen such 
that V x Q has the appropriate singularities along the 
vortex cores, where / vanishes. By minimizing this F 
with respect to ip, A or /, Q, one obtains the GL equa- 
tions together with the appropriate boundary conditions. 
For the superconducting film considered here, one has to 
add the energy of the magnetic stray field outside the 
film, which makes B continuous at the film surface, see 
below. 

The 3D solution of the GL equations for an infinitely 
large, thick or thin film with periodic lattice of vortices 
perpendicular to the film plane, can be obtained numer- 
ically by a modification of the 2D method developed in 
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FIG. 1: Magnetic field lines (top) and profiles of order param- 
eter \ip\ 2 — uj(x,0,Zi) and magnetic field B z (x,0, Zi) (bottom) 
for a superconductor film calculated from Ginzburg-Landau 
theory for the triangular vortex lattice. Shown is the exam- 
ple b = B/B C 2 = 0.04, n = 1.4, triangular lattice with vortex 
spacing (unit length) xi = 3" 1/4 (2$ /B) 1/2 = 5a;i(B c2 ) w 
10A, film thickness d = 0.8zi w 8A. Top: The left half shows 
the field lines that would apply if the field in the film would 
not change near the surfaces z — ±d/2 marked by dashed 
lines. The right half shows the correct solution. The density 
of the depicted field lines is proportional to the local induc- 
tion |B|, see Appendix A and Fig. 2. Bottom: The solid lines 
show ui and B in the center of the film [z = 0) and the dashed 
lines at the film surfaces. The average induction B is marked 
as a dotted line. 



Refs. ITlllTl This is possible since in any plane z = const 
parallel to the film the solutions for the ideal vortex lat- 
tice are still periodic. Actually this periodicity applies 
even to tilted and arbitrarily curved vortex lines, and to 
anisotropic superconductors, which may be computed by 
a similar method. These more complex problems will be 
considered in future work. 

For the present problem of straight vortices along z one 
may choose a general ansatz for the magnitude of the GL 
function f(x,y,z) — \ip(x, y, z)\ in form of the following 
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3D Fourier series for the smooth function / : 
uj(t) = f 2 = ^a K (l - cosK_Lr ± )cos^ 

K 



(3) 



Here r = (x,y,z), r±_ = (x,y), K = (K x ,K y ,K z ), K ± = 
(K x , K y ). In all sums here and below the term = 
is excluded. For vortex positions R = R„ m = ijnx\ + 
nx2, ny2) the reciprocal lattice vectors are Kj_ = K m „ = 
(2ir I S)(my2, 12x1+171x2) with S = x\y2 = §o/B the unit 
cell area and m = 0, ±1, ±2, . . . , n = 0, ±1, ±2, .... 
The z-component of K is chosen as K z — (2n/d)l with 
I = 0, 1, 2, . . . and d the film thickness. This ansatz 
guarantees that /(R, z) = at the (straight) vortex cores 
and that at the two surfaces of the film z = ±d/2 one 
has d(f 2 )/dz — 0, as it follows from the variation of the 
GL free energy functional (2). If only the term K z = 
is kept, the ansatz (3) reduces to that for the 2D vortex 
lattice in Ref. [ij. Formally, the 3D Fourier series (3) may 
also be expressed as a 2D Fourier series with z dependent 
coefficients a~K ± (z) = J2k 8 k cos^z. 

For the supervelocity Q and magnetic induction B = 
V x Q inside the film we chose the general ansatz 



Q(r) = Q A (r x ) + q(r), 

B(r) = Bz + b(r), (b(r)> = 0, 

b(r) = Vxq(r) . 



(4) 



Here Qa(x, y) is the supervelocity of the Abrikosov B c2 
solution, which satisfies 



VxQ y 



B 



R) 



R 



(5) 



where ^(r^) = 8{x)5{y) is the 2D delta function and $ 
the quantum of flux, <I>o = 2tt/k in reduced units. For- 
mula (5) shows that Qa is the velocity field of a lattice 
of ideal vortex lines but with zero average rotation. Near 
each vortex center one has Qa(i"_l) ~ z x r' ±/(2nr' 2 ) 
and /(r) 2 cx r' 2 with r'j_ = rj_ R- QaOm) ma Y be 
expressed as a slowly converging Fourier series by inte- 
grating (5) using divQ = divQ^ = 0. It is, however, 
more convenient to take Qa from the exact relation 



Qa(f_l) = 



Vloa x z 

2 KLUA 



(6) 



where oja( x i v) = f( x -,y) 2 is the Abrikosov B C 2 solution 
given by a rapidly converging series of type (3) with z- 
independent coefficients 
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-(-i) 



eM-K 2 mn S/(87T)} 



and 



(7) 



4 K, 



for general vortex-lattice symmetry, 

— (— \) v exp(— TKvjyA) {v — m 2 + mn + n 2 ) for the 
triangular lattice. The uja from (7) is normalized 
to (ujA{x,y)) = 1; this yields the strange relation 



4 K 



= 1 for any lattice symmetry. Another 



(Vwi/w^) 2 — \7 2 lua/ua — Att/S — const, although both 
terms diverge at the vortex positions; this relation fol- 
lows from (5) and (6) using B = $ /S = 2ir/(nS). The 
useful formula (6) may be proven via the complex B C 2 
solution iPa[x, y); it means that near B C 2 the second and 
third terms in the F, Eq. (2), are equal. 

The general ansatz for q(r) = {q x ,<ly,<lz) is a Fourier 
series for all three components, satisfying V x q = b. 
For simplicity here I shall assume q z — 0, which means 
planar supercurrents. In the considered case of vortices 
perpendicular to the film plane this is an excellent ap- 
proximation, which is exact in the limit of small induc- 
tion and probably also at large inductions B w B c2 , and 
it is exact for thin films. I further assume V • Q = 0, 
which is exact in several special cases (e.g. for B <C B C 2 
and B w B C 2) and is possibly exact even in the general 
case, though I did not find a proof for this. Note also 
that within the circular cell approximation^^ both as- 
sumptions are satisifed. With these two assumptions q 
is fully determined by the z-component of b = (bj_, b z ): 



K 



b z (r) — 5k cos Kj^r^ cos K z z , 
Mr) 
q(r) 



\ ; K ± K Z 

/ ok — TTn — smKir sm K?z 



K 



K 2 



y2 &K Z y2 X sin K -L r X cos K z z , (8) 
Jy i 



K 



with K± = |Kj_|. Formally, these 3D Fourier series (8) 
may also be expressed as 2D Fourier series with z de- 
pendent coefficients &k ± ( 2 ) = J2k bjt cos K z z and their 
derivatives b' K± (z). The solution is thus completely de- 
termined by the two infinite sets of scalar Fourier coef- 
ficients ok and 6k, which are obtained by minimizing 
the total free energy with respect to these coefficients for 
given parameters n and B and film thickness d. For the 
computation I shall use a large but finite number of <zk 
and 6k in the sense of a Ritz variational method. 

The total free energy Ftot per unit volume of the infi- 
nite film is the sum of the GL free energy, Eq. (2), and 
the stray- field energy F stray . In reduced units and re- 
ferred to the state where i\> = and B(r) = Bz = fioH a z 
one has with g = (V/) 2 /k 2 = (Vw) 2 /(4k 2 cj): 



F tot = (-co + V + g + uQ 2 + b 2 ) + 



F stv&y = 2 (B(r) 2 -B 2 ) x , y dz. 

Jd/2 



(9) 



strange property of the Abrikosov solution (7) is that 



The factor of 2 comes from the two half spaces above 
and below the film, which contribute equally to i^tray 
The stray field B(x, y, z > d/2) with constant planar av- 
erage (B(x,y,z)) XtV = Bz is determined by the Laplace 
equation V 2 B = (since V • B = and V x B = in 
vacuum) and by its perpendicular component at the film 
surface z = d/2, since B z has to be continuous across 
this surface. This yields with (8) the stray field: 

B z (x,y,z > d/2) = 



4 



E fe k x cos(Kxrx) exp[-Kx(z - d/2)] + B . 



B ± ( Xl y,z>d/2) = 



K 



E sin(K ± r x ) e X p[-K±(z - d/2)} 



K 



& Ki = *>kA* = d/2) = E b * cos ^) ' ( 10 ) 
i 

(7 = 0, 1, 2, . . .). For spatial averaging we shall need the 
ortho-normality relations valid for =/= 0: 



[ cos(K ± r ± ) cos(K' i r ± ) ) x>v 



;sin(Kxrx)sin(K / _ L rx)) :c ,!/ = o<5klK' x > (11) 



2wlz 2-kI'z \ „ 1 + 5i o 



cos — COS - 

d d 



sin sm 

d d 



= 5u> 



= Su> 



2 ' 
1 - 5, 



(12) 



Averaging the squared stray field over x and y and using 
(11), (12) one obtains terms (6 K J 2 exp[-2K j_(z - d/2)}, 
and thus F s t ray in Eq. (9) becomes 



stray 



E 

K 



'bkj 2 
K ± 



(13) 



The Fourier coefficients ok and 6k may be computed 
by iterating appropriate GL equations as shown in Ref. 
ITllll2l Minimizing F, Eq. (2), with respect to / and Q 
yields the two GL equations for bulk superconductors 



2 V 2 / = -/ + f + fQ 2 



j=VxB=VxVxQ 



-fQ- 



(14) 
(15) 



The first GL equation (14) applies also to our film; the 
second GL equation (15), too, but when it is written as 
an equation for the 6k it has to be supplemented by a 
stray-field term ~ SF s t ia . y /SCl on its r.h.s., which orig- 
inates from the boundary condition for B. A possible 
iteration equation for the ok is obtained from (14) using 
the relation 2/V 2 / = V 2 uj - (X7 lo) 2 / (2lo) to give 



V 2 to = 2k 2 ( — lo + uj 2 + ujQ 2 + g ) 



(16) 



with g = (Vw) 2 /(4k 2 w) as above. Note that V here 
means the 3D Nabla operator, while the similar Eq. (9) 
of Ref.[Tl] is 2D. To obtain better convergence of the iter- 
ation I subtract a term 2k 2 u> on both sides of (16), such 
that K 2 is replaced by K 2 + 2 k 2 ; this choice yields fastest 
convergence. Using the ansatz (3) and the orthonormal- 
ities (11), (12) we then obtain an iteration equation for 
the ok: 



((uj 2 -2lu + ujQ 2 +g) cosK^r^ cos K z z) 
\ (Sk z: o + 1) (K 2 /2k 2 + 1) 



, (17) 



where (...) averages over x,y,z. In particular, if lo and 
Q do not depend on z, Eq. (17) reduces to Eq. (11) of 
Ref. and yields ok = for all K z ^ 0. Other forms 
of iteration equations for the ok are possible, e.g. one 
containing in the denominator K\ instead of K 2 , but 
one should choose that which yields fastest convergence 
of the iteration. In general, the iteration of any equation 
for some parameter a given in the original form a := 
F(a) may be made more stable or faster converging by 
rewriting it in the form a := (1 — c)a + cF(a) with some 
constant c < 1 (or even c > 1 in some cases). 

The convergence is accelerated by alternating the 
iteration step (17) with an iteration step that changes 
only the amplitude of to but not its shape. Namely, from 
dF tot /d{uj) = we obtain 



9 - uQ 2 ) 
(lo 2 ) 



(18) 



Similarly, an iteration equation for the 6k is obtained 
from the equation dF to t /36k = by reordering the terms 
appropriately. From Eq. (8)-(13) one has 



d(b 2 ) = b Jl+5 lt0 )Kl+jl-5 ltQ )K* (jg) 



96 K 

d(coQ 2 ) _ 2P K 



2K\ 



db 



K 2 



1 OF 



K 

.tray _ , s 2 COs(7rZ) 



(20) 
(21) 



d db K dKj_ 
with 6 Ki from Eq. (10) and 

Pk = (io(QyK x — Q x K y ) sin~Kj_r± cos K z z) . (22) 

Equating the sum of the terms (19)-(21) to zero and 
adding and subtracting an appropriate term c(w)6k that 
improves the convergence (with some constant c ~ 1 or 
larger), one obtains an iteration equation for the 6k : 



2P K + c(lo) b 



K 



fiT_L& Ki cos(7r0 



S lfi Kl + i (1 - S h0 )K 2 + c(uj) 



(23) 



The solutions w(r), B(r), and Q(r) are obtained iter- 
atively by first finding the 2D solution as in Ref. EUl 
keeping only the terms with K z = and starting, e.g., 
with ok = (1 — B/B C 2)a^_ ± and 6k = and then it- 
erating the three equations (17), (18), (23) by turns a 
few times; after this, the 3D solution is obtained by con- 
tinuing this iteration with the terms for all K z included 
until the coefficients ok and 6k do not change any more. 
With the empirical choice c ~ 3 + (0.4 + 606 2 )K 2 a;i/d this 
iteration is stable for all 6, k, and d and the free energy 
decreases smoothly until it becomes stationary (with ac- 
curacy 10~ 14 ) after 25 ... 50 iteration steps. 
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FIG. 2: Comparison of the magnetic field lines plotted either 
as stream lines (solid lines) that flow along the exact direc- 
tion of the magnetic field but do not show the correct 2D 
flux density, or as contour lines (dashed lines, see Appendix 
A) that show the correct flux density but have only approx- 
imately the orientation the magnetic field. Shown are the 
examples b = 0.04, triangular lattice, with (right plot) k = 2, 
d = 0.6xi « 4A and (left plot) « = 1, d = 0.6xi « 8A. The 
horizontal dashed lines indicate the film surfaces z = ±d/2. 
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FIG. 3: The magnetic field lines of the vortex lattice in films 
of various thicknesses d/X w 4, 2, 1, 0.5, corresponding to 
d/xi = 0.4, 0.2, 0.1, 0.05, for b = 0.04 and k = 1.4. Depicted 
is the field in the plane y = in one lattice cell. The film 
surfaces are marked by two dashed lines. The field lines of an 
isolated vortex in such films are shown in Fig. 2 of Ref. 9. 



FIG. 4: The variances of the longitudinal and transverse 
components of the magnetic induction, a z (z) and <r±(z), de- 
fined by Eqs. (24), (25), plotted versus z/xi (xi = vortex 
spacing) at reduced induction b — 0.04 for films of various 
thicknesses d/xi = 0.1, 0.2, 0.4, 0.6, 0.8, 1, 1.2, and 1.4. Top: 
For k = 1, yielding X/xi = 0.0743. Bottom: For k = 3, thus 
X/xi = 0.223. While <j z (z) decreases monotonically with in- 
creasing 2, u±(z) has a sharp peak at the film surface z — d/2. 
Outside the film (\z\ > d/2) one has a z (z) = a±(z). 



III. SOME RESULTS 
A. Magnetic field and order parameter 

Figure 1 shows one example for the resulting mag- 
netic field lines and some cross sections of lo(x, y, z) and 
~B(x,y, z) along x in the plane y = at z — (center 
plane of the film) and z = d/2 (film surface), for a film 
of finite thickness d = 0.8x\ rts 8A at reduced induction 
b = B/B C 2 = 0.04 and GL parameter k = 1.4, yield- 
ing for the triangular vortex lattice a vortex spacing of 
x x (B) = 5xi(S c2 ) = 1.25d w 10A. The left half of Fig. 1 
shows the field lines that result if the unchanged 2D bulk 
solutions for B{x,y) and cj(x,y) are assumed inside the 
film. The right half shows the correct solution, exhibit- 
ing smooth field lines across the film surface, and a very 
weak widening of the vortex core near the surface. 

Figure 2 shows the magnetic field lines for a film with 
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FIG. 5: As Fig. 4 but for k, = 1.4 and two b values. Top: 
Large b = 0.3, yielding X/xi — 0.285. Bottom: Low b = 0.01, 
thus A/asi = 0.0520. 



thickness = 0.6xi at b = 0.04 for k = 2 (rf w 4A, 
left) and k = 1 (d w 8A, right), triangular lattice. The 
solid lines are the stream lines of B(x, 0, z) = 0, B z ); 
they have the correct slope of B and start at equidistant 
points far away from the film surface, where B = B Q = 
(0, 0, B) = constant (in infinitely extended films the aver- 
age induction B equals the applied induction B a outside 
and inside the film), but their ID density is not propor- 
tional to the 2D flux density B = |B|. The dashed lines 
in Fig. 2 are field lines that have approximately the slope 
of B(x,0, z) and have a density proportional to B, see 
Appendix A. This type of field lines is depicted also in 
Figs. 1 and 3. 

In Fig. 3 the magnetic field lines are shown for films of 
various thicknesses d/xi = 0.4, 0.2, 0.1, 0.05 for b = 0.04 
and k — 1.4 as in Fig. 1, where djx\ = 0.8. These 
thicknesses correspond to d/Aw 4, 2, 1, 0.5 (and 0.25 in 
Fig. 1). At low inductions b <C kT 2 and not too small 
K > 5, these field patterns may also be obtained by lin- 
ear superposition of the fields of isolated London vortices 
given by Eqs. (5)-(9) of Ref.0, with appropriately cut-off 
vortex core introduced to consider the finite coherence 
length £. This superposition also applies to nonperiodic 
vortex arrangements. 
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FIG. 6: The spatial variance (Abrikosov parameter) 
[3(z) — { uj 2 }x, v /{ u }x,yi °f the order parameter u(x,y, z) = 
\tp(x,y,z)\ 2 , plotted as W 4 [f3{z) - /3(0]] + 100[/3(0) - 1] ver- 
sus z/(d/2) for several reduced inductions b — 0.99, 0.8, 0.4, 
0.2, 0.1, 0.05, 0.02, and 0.01 for k = 1.4 and for both trian- 
gular (solid lines) and square (dashed lines) vortex lattices. 
Top: Thick film with d — 4x C 2. Bottom: Thinner film with 
d — x C 2- Here x C 2 = x\\/b w 2.7^ is the vortex spacing 
at B = B C 2- The extremely small variation of f3(z) is en- 
larged by plotting 10 4 x [f3(z) — /3(0)]. Adding the constant 
10 2 [/3(0) - 1] allows to identify (at z = 0) the bulk Abrikosov 
values /3(0) w /3 A = 1.1596 (1.8034) occurring for the trian- 
gular (square) lattice in thick films when b 1. 

B. Variance of the magnetic induction 

Figures 4 and 5 show the two relative variances a z and 
<t_l of the magnetic induction defined here as 

a z (z) = ([B z (x,y,z)-B} 2 )l^/B (24) 
a ± (z) = (B x (x, y ,z) 2 +B y (x 1 y,z) 2 ) 1 J*/B. (25) 

These measures of the relative variation of the periodic 
induction depend on z: Deep inside thick films, a z (z) 
reduces to the variance of the 2D vortex lattice in the 
bulk, (Tbuik, computed, e.g., in Ref. as function of b 
and k, and one has a±_ = since B||z for the considered 
case. As one approaches the surface from inside the film, 
a z (z) decreases and erj_(z) increases until they coincide 
at the surface z — d/2. Outside the superconductor one 
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FIG. 7: The energy Fstray of the magnetic stray field (solid 
lines) and the surface energy F sur f of the film (dashed lines) 
plotted versus the film thickness d for k = 2 and b = 0.02, 
0.15, and 0.4, see text. These energies per unit area are plot- 
ted in units /io-ffcC/20000 to enlarge them to order of unity 
and show their close coincidence at large b. 



exponentially with z, a z = cr± oc exp(— .?Tio|z|), where 
if io w 2-k/xi is the shortest reciprocal lattice vector of 
the vortex lattice defined below Eq. (3). Thus, <r z (z) 
decreases monotonically from its bulk value Obuik inside 
a thick film to zero far away from the film, reaching at 
the surface approximately half its bulk value (for thick 
films). In contrast to this, the transverse variance cr±(z) 
increases from zero at z — and reaches a sharp cusp- 
shaped maximum at the surface, where it joins cr z (z) and 
then decreases again to zero away from the film. For large 
k > 2 and d/x\ > 0.7 these curves are approximately 
symmetric, 

<t z (\z\ < d/2) = (Tbuik - <r z (d- \z\) 
a±(\z\<d/2)na z (d-\z\) 
a±(d/2) = a z (d/2) w ^bulk/2 . (27) 

This is so since for A> xx/2-k the outer and inner mag- 
netic stray fields are symmetric. 

Shown in Figs. (4) and (5) are both variances for var- 
ious values of b and k for 8 film thicknesses d/x\ = 0.1, 
0.2, 0.4, 0.6, 0.8, 1, 1.2, and 1.4 for the triangular vortex 
lattice (solid lines) and for d/x\ = 0.1, 0.4, 0.8, and 1.2 
for the square vortex lattice (dashed lines) . The variances 
for these two lattice symmetries are almost identical. 
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FIG. 8: The thick-film limits of the stray-field energy F stra y 
(solid lines) and the surface energy F sur f (dashed lines) plotted 
versus the reduced induction b for ft = 0.71, 2, and 20 in units 
MoHcC/2- 10 4 . 



has exactly 

<j 2 z (z > d/2) = a\{z > d/2) 
= E( fe K J 2 exp[-2if ± (z - d/2)} . (26) 



K 



This follows from Eqs. (10) for the magnetic stray field. 
At large \z\ — d/2 » xi/(2w), the variance decreases 



C. Variance of the order parameter 

Figure 6 shows the variance of the order parameter 
u>(x, y, z) — \ip(x, y, z)\ 2 inside films with periodic vortex 
lattice, 



(<jj(x,y,z) 2 ) Xt y 
{u(x,y,z)) 2 ' 



(28) 



At large reduced inductions b — B /B C 2 — > 1 in the middle 
of thick films (3(z) coincides with the Abrikosov param- 
eter Pa = 1.15960 for the triangular vortex lattice and 
[3a — 1.18034 for the square lattice. At low inductions 
6^0 one has /? — > 1 since the order parameter uj is con- 
stant except in the small vortex cores. Figure 6 shows 
that P(z), and thus the shape of ui{x,y, z), changes very 
little with z. For films of thickness d > 4x C 2, one has 
constant bulk [3 in a central region around z = 0, and 
as z approaches the surface of the film, (3(z) increases 
by at most a factor of 1.0033 within a layer of thickness 

x C 2- Here x C 2 = x\{B C 2) = x\\fb is the vortex distance 
at B — B C 2- One has x 2 c2 = (47r/v / 3)^ 2 for the triangu- 
lar and x 2 2 = 2tt^ 2 for the square vortex lattice. The 
maximum change occurs at b w 0.15 and d > Ax C 2- For 
thinner films and larger and smaller 6, the variation of 
(3(z) is even smaller. 

Thus, to a very good approximation, one may assume 
that the order parameter u)(x, y, z) of the vortex lattice 
inside films in perpendicular magnetic field is indepen- 
dent of z and for not too thin films has the same form 
as the 2D order parameter u>(x, y) of the bulk vortex lat- 
tice. For very thin films with d <C £ at b <C 1 the vortex 
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cores are slightly wider than in the bulk. For example, 
at b = 0.04, k = 1.4, the core width increases by about 
25 % when d/£ decreases from 0.5 to 0.005, but then 
saturates and does not increase further in thinner films. 
This is just the interval of d in which the modulation 
1 — B mul / B max of the periodic magnetic field B(x,y,0) 
decreases from ss 1 to -C 1 since the effective penetration 
depth 2X 2 /d becomes larger. 

Close to B C 2 the constancy along z of the GL function 
ip(x,y,z) applies to thicker and thicker films. This nu- 
merical result is consistent with the finding in Ref. .24 of 
a correlation length l z = £/(2\/l — b) that diverges for 
6 — > 1 and describes the extension along the vortex lines 
of perturbations in ip(x, y, z) caused by small material 
inhomogeneities (pins). Interestingly, a similar diverging 
length ^/(2-vT — b) describes the long axis (along z) of a 
cigar-shaped superconducting region (nucleus) that nu- 
cleates at applied fields above B C 2 at a small defect with 
transition temperature T c (r) higher than the bulk T c i& 
In superconducting films of thickness d < 1| -1 / 2 , or 
at applied fields satisfying \B a /B C 2 — 1| < £, 2 /d 2 , small 
inclusions or precipitates are thus expected to cause a 
virtually 2D perturbation that has no z-dependence. 



D. Surface energy 

Next I consider the additional free energy caused by 
the presence of the two surfaces of the film. This energy 
per unit area of the film is composed of two terms: 

(a) i^tray , the magnetic energy of the stray field out- 
side the film, defined by Eq. (9) and expressed in Eq. (13) 
in terms of the Fourier coefficients bfj of the field com- 
ponent B z (x,y,d/2) at the surface, 

(b) F sur f, the actual surface energy defined as the dif- 
ference of the free energy of the film per unit area, F 3 £,d, 
minus the 2D bulk free energy density of the infinite vor- 
tex lattice, i*2D, times d, thus 

Furf = (F 3D ~ F 2D )d. (29) 

The total surface energy, originating from both sur- 
faces of the film, is the sum of these two terms, F s t r ay + 
-Fsurf • Both terms tend to a constant when the film thick- 
ness d increases above the vortex spacing x\ . These thick- 
film values of -F stra y and F sur f are of the same order, and 
they are approximately equal for large k> 1 and also at 
large reduced inductions 6 — > 1. This is so since the or- 
der parameter lo(x, y, z) in the film is nearly independent 
of z, and thus F sur f is virtually only of magnetic origin, 
i.e., it is the energy of the magnetic field change caused 
inside the film by the surface. When the magnetic pene- 
tration depth is large, A > xi/2tt <C d, this "inner stray 
field" is symmetric to the outer stray field. This equality 
applies also near 6=1, since inside the superconductor 
the magnetic screening is reduced by the reduction of the 
order parameter and thus the effective penetration depth 
A' = A/» 1/2 ~ A /(l - b) 1 ' 2 increases-^ 



The dependence of -F s tray and F sur f on the film thick- 
ness is shown in Fig. 7 for k = 2 and b = 0.02, 0.15, 
and 0.4. With increasing d both energies increase from 
zero and saturate to constant values at about d/x\ > 0.7 
for all k and b. Figure 8 shows these thick-film limits of 
-Fstray and -Fsurf as functions of the reduced induction b. 
Note that F stra y is slightly larger than F sur f ; this reflects 
the fact that the stray field inside the film is screened 
by A' < oo, while outside the film there is no screening 
(A = oo). Both F s t r ay and F sur [ vanish at b — > and 
at b — > 1 and have a maximum in between. At b — > 
one has Ftray ~ F sm f oc b since each vortex contributes 
separately. At b — > 1 one has F stra y ~ F snI { oc (1 — b) 2 / n 
(in units ^,qH 2 X) since the amplitude of the periodic 
B z (x,y,d/2) decreases as 6^ oc (1 — b)/n and the depth 
of the stray field is Kio w xi/2ir oc £ = X/k. Therefore, 
when plotted in units /io-ff 2 £, all curves i^tray and all 
F sur f practically coincide for all n > 1 and b > 0.4, see 
Fig. 8. 

Note that the total surface energy F stlay + F snI { is very 
small, never exceeding the value 8 • 10~ 4 /io-ff 2 £ reached 
at k> 1 and b w 0.2/ y/K. 



E. Shear modulus of vortex lattice 

Finally, the elastic shear modulus cqq of the vortex 
lattice will be discussed. As shown in Ref. 0, the shear 
modulus of the triangular vortex lattice can be expressed 
with high accuracy by the difference of the free energies 
of the rectangular lattice, F rec t (with x 2 — and y 2 = 
y/3xi/2), and the triangular lattice, F tI (with x 2 — X\j2 
and same y 2 = V3xi/2), 

c 66 = (3^ 2 /2)(F.ect-F tr ). (30) 

This is so since the free energy for constant unit cell 
height ?/2 varies practically sinusoidally with x 2 : F(x 2 ) ~ 
Ftr + [ l+cos(27ra;2/xi) ](F rect — F tr )/2, thus the definition 
c 66 = d 2 F/da 2 at small shear angle a = arctan[(a;2 — 
0.5) jy 2 \ yields Eq. (30)i^ Expressed in units HoH 2 , 
the shear modulus depends on three variables: cqq — 
cee(b, k, d). There are thus several ways to present the 
numerical data for C66, each yielding different physical 
insight. 

One result is that in the limit of small film thickness 
d <C £ the shear modulus for a film with any k tends to 
the bulk shear modulus at k — > oo, as already obtained 
by Conen and Schmid. 14 This finding may be understood 
from the fact that in thin films the vortices are Pearl vor- 
tices that have a long interaction range 2 A 2 jd exceeding 
the London penetration depth Xw^ This argument yields 
the correct limit n e g — > oo for d/X — > 0, but for k < oo 
the C66 of films does not quantitatively coincide with the 
bulk C66 for an effective k d s = 2X 2 /d£ — (2X/d)n ^> k, 
since cqq is determined not only by the range but by the 
full shape of the interaction potential between vortices, 
which differs for Abrikosov 1 and Pearl 4 vortices. If this 
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FIG. 10: The shear modulus cg6 of the triangular vortex 
lattice in films of various thicknesses like in Fig. 9 bottom, 
but for k = 1 (top) and k = 2 (bottom), in units noHc /1000. 



FIG. 9: Top: The shear modulus cg6 of the bulk (d —> oo) 
triangular vortex lattice as function of the reduced induction 
b = B/B c2 for GL parameters k = 0.4, 0.5, 0.6, .707, 0.75, 1, 
1.4, 2, 3, 5, 7, 10, 100, in units fi H^/1000. For k < 2" 1/2 = 
0.707 one has formally C66 < 0, though vortices and a vortex 
lattice are unstable in such type-I superconductors. Bottom: 
The shear modulus C66 of the triangular vortex lattice in films 
with thicknesses d/£ = 0.1, 0.32, 0.56, 1, 1.8, 3.2, 5.6, 10, and 
32, plotted versus b for k = 0.5. This C66 is positive, i.e. the 
triangular vortex lattice is stable, for sufficiently thin films or 
small inductions. For d > ( the bulk ces at the same k = 0.5 
is reached (dash-dotted line), and for d <C £ the bulk cqq in 
the limit k> 1 (dashed line) . 

potential is V(r) with r = (x 2 + y 2 ) 1 ^ 2 and the vortex 
density is n v = B/&q, one has for a 2D triangular lattice 
with positions R mn defined below Eq. (3)i2£i2£ 

c 66 = Jj. [ y" {Rmn ) + 3R mn V' (R mn ) ] , (31) 
m.n 

see also Eqs. (9) and (11) of Ref.El 

Figures 9 and 10 show cqq as a function of the reduced 
induction b for different film thicknesses d expressed in 



units of the GL coherence length £, d/£ = lO^" 4 —' 6 )/ 4 = 
0.1, . . ., 32, and for GL parameters k = 0.5, 1, and 2. For 
large k > 5, the curves C6e(b) for various thicknesses are 
very close to each other and to the bulk cqq. In general, 
the curves for finite film thickness all fall between the two 
limiting cases d — > oo coinciding with the bulk cqq(k), 
and d — ► coinciding with the bulk Cqq(k — oo). This 
interval is very small for large k and not too small b since 
C66( K = oo) ~ cqq(k) oc kT 2 . This means that for large 
k > 5 the shear modulus is nearly the same for thin and 
thick films. Note that for the bulk and k > 5 one has 
c 66 w B<S> Q /(16TTfi \ 2 ) oc b for 1/(2k 2 ) < b < 0.15 (see 
Fig. 9 top); this applies also to films. For n < 5 we 
confirm the finding of Refi^ that cee oc 6 3 / 2 for b <C 1, 
but this law applies only to intermediate film thicknesses 
0.5 < d/X < 3 at b < 0.1. 

An interesting feature can be seen from Fig. 9. The 
upper part shows the bulk Cee(b,K) for values k = 0.4 
to oo, i.e., also for k < l/y/2 — 0.707 corresponding 
to type-I superconductors, in which the vortex lattice is 
energetically unfavorable in the bulk. For k < 0.707 one 
finds negative cqq < 0. This means the bulk rectangular 
and square vortex lattices^ have lower energy than the 
triangular lattice (the Meissner state without vortices has 
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FIG. 11: Dependence of the shear modulus egg of the trian- 
gular lattice on the film thickness d. Plotted is the additional 
rigidity caused by the film surfaces in form of f{d, b, k) = 
[cee(d,b, k) — C6@(oo,b, «)] ■ (d/£) ■ k versus d at b = 0.3 for 
various k = 1 ... 20 (solid lines, cee in units fioH 2 /1000). 
The dashed lines show the fit, Eq. (32), with I = £. 



even lower energy in this case). However, as can be seen 
in the lower plot for films with n — 0.5, in sufficiently 
thin type-I superconductor films the triangular vortex 
lattice can be stable (i.e., cqq > 0) when the induction is 
sufficiently small. This behavior was seen also in Ref. 

The dependence of cqq on the film thickness d is visu- 
alized in Fig. 11 by plotting the difference 

f(d, b, k) = [c 66 (d, b, k) - c 66 (oo, b, k) ] ■ (d/£) ■ k (32) 

(an energy per unit area) versus d at various k values for 
b = 0.3 (near the maximum of cqq). One can see that 
this function saturates when the film thickness exceeds a 
few coherence lengths £. For all values of k and b > 1/k 2 
one can fit these curves closely by 



/(d)w/(oo)[l-exp(-d/0; 



(33) 



In Fig. 11 (at b — 0.3) the length I of the best fit acciden- 
tally coincides with £, but in general I depends on b and 
is proportional to the vortex spacing x\. I ~ 0.195:ci = 
(v%/4tt)^i = V2/K 10 , thus l/Z sa 3 1 / 4 (2tt)- 1 /2/^5 = 
0.525/V^, yielding I = 0.96^ at b = 0.3. This saturation 
means that the additional rigidity of the vortex lattice 
caused by the film surfaces and measured by f(d,b,K), 
becomes independent of d in films thicker than a few co- 
herence lengths, and thus one has cee(d) — C66(oo) oc 1/d. 
For thin films with d <C £ one has f(d, b, k) oc d since 
c 66(^) — C66(oo) is a constant independent of d. 

Figure 12 visualizes the k dependence of the 
shear modulus of the triangular lattice by plotting 
[c66(oo, 6, oo) — cee(d,b, k)] ■ k 2 versus 1/k at b = 0.1, 
0.3, and 0.7, for film thicknesses <i/£ = 0.1 ... 10. Also 
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FIG. 12: Dependence of the shear modulus c^a of the trian- 
gular lattice on the GL parameter k. Plotted is the function 
[c66(oo, b, oo) — C(i6(d,b, jc)] ■ k 2 versus 1/k at b — 0.1 (top), 
0.3 (middle), and 0.7 (bottom), for various film thicknesses 
d/£ (c66 in units fioH 2 /1000). The dashed lines show the 
limits d = (z-axis) and d = oo (upper line) coinciding with 
[cg6 lk («; = oo) - C6e lk (K)] • k 2 . It is clearly seen that the dif- 
ferences of any two ce>6 values vanish as k~ 2 when k — > oo. 
This asymptotic law is good even for k > 2, and it practically 
applies to all k > 0.71 at large inductions b > 0.7. 



shown are the limits d = (the x-axis), proving that 
Cm{d — > 0, 6, k) = cee(d — > oo,6, k — > oo) for any k, 
and d = oo (upper dashed line) that coincides with 



„bulk 



(« = «>)- c 6 b 6 ul >) 



. These plots prove that 
the differences of any two cqq values vanish as k~ 2 when 
k — ► oo. One can see that this asymptotic law is a good 
approximation even for not so large k > 2, and at large 
b > 0.7 it practically applies to all k > 0.71, and even for 
type-I superconductors with not too small k < 0.71. 

In Fig. 13 the numerical cqq{cI, b, k) is compared with 
the analytical result of Conen and Schmid, Fig. 1 of 
Ref. 0, valid at large inductions 1 — b <C 1 . Their re- 
sult was obtained from an elegant expression derived by 
Lasher— for the free energy of films with vortex lattice of 
arbitrary symmetry at b — > 1. Lasher 1 - implicitly used 
the fact that the magnetic stray field inside the film is 
not screened in this limit of b — * 1 . 
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FIG. 13: The shear modulus cse oc (1 — b) 2 of the vortex 
lattice in films at high inductions 1 — b <C 1, plotted versus 
the film thickness d. For comparison with Fig. 1 of R,ef.H3. C66 
is in the same units (1 — b) 2 fioH 2 /2. The horizontal dashed 
lines denote the bulk C66- This plot, computed for b = 0.98, 
applies approximately too all b > 0.6 and ag rees with the 
analytic expression plotted in Fig. 1 of Ref. UA . 
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FIG. 14: The magnetization of infinite films with thickness 
d/£ = 0.1, 1 ,3, 10, oo, containing a triangular vortex lattice 
with one flux quantum per vortex. Plotted is — m — —M/H^ 
versus h = H/H c2 = b = B/B c2 for k = 0.5, 0.707, 1, 1.5. 



F. Magnetization of infinite films 

For infinitely extended films the average magnetic in- 
duction B equals the applied field, B = ^H a , and the 
magnetization M is defined as the magnetic moment 
per unit volume of the film. The demagnetization fac- 
tor of infinite films is A = 1, and thus the effective 



lower critical field at which the first vortices penetrate is 
H' cl = (1 — N)H C \ = 0. For the magnetization of super- 
conductors wit h ge neral demagnetizing factor < N < 1 
see, e.g., Refs. [1912a Noting that the total free energy 
per volume F tot , Eq. (9), equals the work done by the 
applied field, which may be written as — J M dB, one 
obtains the relation M — —dF tot /dB. Figure 14 shows 
magnetization curves for films of various thicknesses d/£ 
= 0.1, 1, 3, 10, and oo for GL parameters k = 0.5, 
1, and 1.5. 

For thick films (d ^> £, A) at k — 1/V2 one has 
Ftat — "5(1 — fr) 2 an d thus —to = l — b= 1 — h; for larger 
k > 0.707 the thick film —m(h) is concave (has positive 
curvature); and for smaller n < 0.707 (type I supercon- 
ductors) —m(h) is convex (has negative curvature) and 
the vortex lattice is not stable. However, for sufficiently 
thin films, even when k < 0.707 the curvature of —m(h) 
can be positive and even may change sign at a certain 
value of h = b. Note that the plotted curves —m(h) for 
various d/£ cut each other approximately at h = b ~ 1/k. 
For thick films the initial slope is —rn'(h)\h = Q = —1 for 
all k, and — m(0) = h cl — H c i/H c2 . This is so since when 
surface contributions may be disregarded, one has for su- 
perconducting ellipsoids — M = H c \ at H = H' cl where 
5 = 0. The enhancement of — m(h) at small h < n for 
thin films in Fig. 14, originates from the energy of the 
magnetic stray field, which enhances the self energy of a 
vortex beyond its bulk value d^oH c \. More details about 
this will be published elsewhere. 



IV. SUMMARY 

It is shown how the Ginzburg-Landau equations can be 
solved for a periodic lattice of Abrikosov vortices in su- 
perconducting films in a perpendicular magnetic field. As 
illustration how well this iteration method works, some 
results are presented. The widening of the magnetic field 
lines as they exit the film surface is correctly obtained, 
Figs. 1, 2, but this leads only to a very small correction 
of the order parameter near the surface, Fig. 1. The vari- 
ance of the transverse component of the magnetic induc- 
tion is sharply peaked at the surface and vanishes deep 
inside and far outside the film, Figs. 4, 5. The variance 
of the periodic order parameter (Abrikosov parameter (3) 
varies very little across the film thickness, by at most a 
factor 1.0033, Fig. 6. The surface energy saturates for 
large film thickness d and vanishes linearly at small d, 
Fig. 7. For not too thin films the surface energy origi- 
nates mainly from the magnetic stray field, which comes 
in approximately equal parts from outside and inside the 
film, in particular for large k or large b, Fig. 8. For 
very thin films the stray field energy may be disregarded 
and the very small surface energy comes mainly from the 
small modification of the order parameter, Fig. 7. The 
shear modulus CQ§{d, b, k) of the triangular vortex lattice 
in thin films approaches the C66(oo,6, 00) of thick films 
(bulk limit) at k — * 00, Fig. 10, the difference being pro- 



12 



portional to k~ 2 , Fig. 12. While the bulk c^(oo : b, n) is 
negative in type-I superconductors (k < 0.707), the cqq 
of sufficiently thin films can be positive and may change 
sign at some value of 6, Fig. 9. The magnetization curves 
of films with k < 0.707 may have positive or negative cur- 
vature, depending on the film thickness, Fig. 14. More 
results will be published elsewhere. Extensions of this 
method to vortices with several flux quanta and to the pe- 
riodic lattice of curved vortices in superconducting films 
in a tilted magnetic field are underway. 



APPENDIX A: PRESENTATION OF FIELD 
LINES 

A practical question is how to plot the magnetic field 
lines of this 3D problem such that they have the correct 
orientation of B(x, y, z) and their ID density (reciprocal 
distance) in the plotted plane is proportional to the mag- 
nitude |B(:r, y,z)\. A simple consideration shows that 
this is possible only for 2D planar problems, when the 
field lines coincide with the contour lines of the vector 
potential, e.g., A y (x, 0, z). But for 3D magnetic fields, 
and even for cylindrically symmetric fields, such 2D plots 
of the field lines are not possible since the magnitude 
|B(x, y, z)\ here is proportional to the 2D density of the 
3D field lines, but not to the ID density of the plotted 2D 



field lines. For our 3D problem of a thick film with a 2D 
periodic vortex lattice we have two possibilities to plot 
field lines that approximately have the above mentioned 
properties. 

First method: One may use numerical programs that 
plot the field lines (stream lines) of the 2D planar field 
B(x,0, z) — (B Xl B z ) [or any other planar cross sec- 
tion of B(ar, y, z)] starting from equidistant points [x — 
Xi,y = 0, z = —zq) far away from the film surface so 
that B(a;, y, z) ~ Bz = const. Such field lines have the 
correct slope, but their density is only approximately pro- 
portional to |B(ar, y, z)\. 

Second method: In this paper the 2D plots of the mag- 
netic field lines show the contour lines of the function 

(p(x,z)= / B z (x, 0, z) dx I j B z (x,0,z)dx, (Al) 

which ranges from ip(0, z) — at x — (vortex cen- 
ter) to (p(x±/2, z) — 1 (middle plane between two vor- 
tices) and has a periodic derivative. Such field lines have 
a ID density along x proportional to \B z {x, 0, z)\, and 
a density perpendicular to these lines which is close to 
|B(a;, y, z)\, since their orientation is close to the orien- 
tation of B(a;, y, z). Figure 2 shows that these two types 
of field lines are very similar. In particular, the contour 
lines of ip(x, z), Eq. (1), have slopes that are close to the 
correct slope. 
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